## This code creates Table 2, Figure 7 (b) and Figure 7 (d)

neighbor_dist_cat <- seq(10,100,10)

# Table 2, 1912 results

data <- colonial_towns_shp[colonial_towns_shp$neighbor_dist <= 50 & colonial_towns_shp$year==1912,]

m_police_specs <- c("m_police_OLS_BAGO_NoGEO","m_police_OLS_BAGO_NoFE","m_police_OLS_BAGO_FULL","m_police_OLS_BAGO_ih","m_police_OLS_BAGO_sd")

m_police_OLS_BAGO_NoGEO <- conleyreg(m_police~replaced+log(houses),data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
m_police_OLS_BAGO_NoFE <- conleyreg(m_police~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
m_police_OLS_BAGO_FULL <- conleyreg(m_police~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
m_police_OLS_BAGO_ih <- conleyreg(m_police_ih~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
m_police_OLS_BAGO_sd <- conleyreg(m_police_sd~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)

m_police_table <- matrix(nrow = 9,ncol=6)

for(i in 1:length(m_police_specs)){
  m_police_table[1,i+1] <- m_police_specs[i]
  model <- get(m_police_specs[i])
  coef <- round(model$coefficients[2],digits=4)
  se <- round(model$coefficients[2,2],digits=4)
  m_police_table[2,i+1] <- coef
  m_police_table[3,i+1] <- paste0("(",se,")")
  m_police_table[4,i+1] <- paste0("[",round(coef-1.96*se,digits=4),",")
  m_police_table[5,i+1] <- paste0(round(coef+1.96*se,digits=4),"]")
  m_police_table[6,i+1] <- nobs(model)
  m_police_table[7,i+1] <- round(model$adj.r.squared,digits=4)
}

xtable(m_police_table)

# Figure 7 (b)

m_police_BAGO <- data.frame(coefficient=double(),se=double(),upper=double(),lower=double())

for(j in 1:length(neighbor_dist_cat)){
  model <- conleyreg(m_police_sd~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,dist_cutoff=100,lat="latitude",lon="longitude",
                     data=colonial_towns_shp@data[colonial_towns_shp@data$year==1912 & colonial_towns_shp@data$neighbor_dist<=neighbor_dist_cat[j],])
  m_police_BAGO[j,"coefficient"] <- model[2,1]
  m_police_BAGO[j,"se"] <- model[2,2]
  m_police_BAGO[j,"upper"] <- m_police_BAGO[j,"coefficient"] + (1.96*m_police_BAGO[j,"se"])
  m_police_BAGO[j,"lower"] <- m_police_BAGO[j,"coefficient"] - (1.96*m_police_BAGO[j,"se"])
  m_police_BAGO[j,"upper_10pc"] <- m_police_BAGO[j,"coefficient"] + (1.645*m_police_BAGO[j,"se"])
  m_police_BAGO[j,"lower_10pc"] <- m_police_BAGO[j,"coefficient"] - (1.645*m_police_BAGO[j,"se"])
  m_police_BAGO[j,"pch"] <- ifelse(sign(m_police_BAGO[j,"upper"])==sign(m_police_BAGO[j,"lower"]),"black",ifelse(sign(m_police_BAGO[j,"upper_10pc"])==sign(m_police_BAGO[j,"lower_10pc"]),"gray","white"))
}

x.axis <- rep(1:10)
labels <- seq(0,100,10)
labels.rev <- rev(labels)
png(file=paste0("Plots/m_police_BAGO_1912.png",sep=""), width=6.5, height=6.5, pointsize=9, units="in", res=600)
par(mar=c(6, 6, 3, 3) + 0.1)
plot(c(0.5,10.5),c(-0.5,0.5),type="n",axes=F,
     xlim=range(0.5:10.5),
     ylim=range(-0.5:0.5),
     xlab="Distance from a location in sittan (km)",
     ylab="Coefficient of headman replacement",cex.lab=1.5,cex.axis=1.5)
box()
axis(2, at=seq(-0.5,0.5,by=0.1), labels=c(-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1))
axis(1, at=seq(0,10), las=2, labels=labels)
segments(x.axis,m_police_BAGO$lower/2,x.axis,m_police_BAGO$upper/2)
points(x.axis,m_police_BAGO$coefficient/2, pch=21, bg=m_police_BAGO$pch)
abline(h=0, lty=3)
dev.off()

library(xtable)

neighbor_dist_cat <- seq(10,100,10)

# Figure 7 (d)

data <- colonial_towns_shp[colonial_towns_shp$neighbor_dist <= 50 & colonial_towns_shp$year==1924,]

data$longitude_2 <- data$longitude

data$latitude_2 <- data$latitude

m_police_specs <- c("m_police_OLS_BAGO_NoGEO","m_police_OLS_BAGO_NoFE","m_police_OLS_BAGO_FULL","m_police_OLS_BAGO_ih","m_police_OLS_BAGO_sd")

m_police_OLS_BAGO_NoGEO <- conleyreg(m_police~replaced+log(houses),data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
m_police_OLS_BAGO_NoFE <- conleyreg(m_police~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
m_police_OLS_BAGO_FULL <- conleyreg(m_police~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
m_police_OLS_BAGO_ih <- conleyreg(m_police_ih~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
m_police_OLS_BAGO_sd <- conleyreg(m_police_sd~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)

m_police_table <- matrix(nrow = 9,ncol=6)

for(i in 1:length(m_police_specs)){
  m_police_table[1,i+1] <- m_police_specs[i]
  model <- get(m_police_specs[i])
  coef <- round(model$coefficients[2],digits=4)
  se <- round(model$coefficients[2,2],digits=4)
  m_police_table[2,i+1] <- coef
  m_police_table[3,i+1] <- paste0("(",se,")")
  m_police_table[4,i+1] <- paste0("[",round(coef-1.96*se,digits=4),",")
  m_police_table[5,i+1] <- paste0(round(coef+1.96*se,digits=4),"]")
  m_police_table[6,i+1] <- nobs(model)
  m_police_table[7,i+1] <- round(model$adj.r.squared,digits=4)
}

xtable(m_police_table)

# Whisker plots showing cutoff points

m_police_BAGO <- data.frame(coefficient=double(),se=double(),upper=double(),lower=double())

for(j in 1:length(neighbor_dist_cat)){
  model <- conleyreg(m_police_sd~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,dist_cutoff=100,lat="latitude",lon="longitude",
                     data=colonial_towns_shp@data[colonial_towns_shp@data$year==1924 & colonial_towns_shp@data$neighbor_dist<=neighbor_dist_cat[j],])
  m_police_BAGO[j,"coefficient"] <- model[2,1]
  m_police_BAGO[j,"se"] <- model[2,2]
  m_police_BAGO[j,"upper"] <- m_police_BAGO[j,"coefficient"] + (1.96*m_police_BAGO[j,"se"])
  m_police_BAGO[j,"lower"] <- m_police_BAGO[j,"coefficient"] - (1.96*m_police_BAGO[j,"se"])
  m_police_BAGO[j,"upper_10pc"] <- m_police_BAGO[j,"coefficient"] + (1.645*m_police_BAGO[j,"se"])
  m_police_BAGO[j,"lower_10pc"] <- m_police_BAGO[j,"coefficient"] - (1.645*m_police_BAGO[j,"se"])
  m_police_BAGO[j,"pch"] <- ifelse(sign(m_police_BAGO[j,"upper"])==sign(m_police_BAGO[j,"lower"]),"black",ifelse(sign(m_police_BAGO[j,"upper_10pc"])==sign(m_police_BAGO[j,"lower_10pc"]),"gray","white"))
}

x.axis <- rep(1:10)
labels <- seq(0,100,10)
labels.rev <- rev(labels)
png(file=paste0("Plots/m_police_BAGO_1924.png",sep=""), width=6.5, height=6.5, pointsize=9, units="in", res=600)
par(mar=c(6, 6, 3, 3) + 0.1)
plot(c(0.5,10.5),c(-0.5,0.5),type="n",axes=F,
     xlim=range(0.5:10.5),
     ylim=range(-0.5:0.5),
     xlab="Distance from a location in sittan (km)",
     ylab="Coefficient of headman replacement",cex.lab=1.5,cex.axis=1.5)
box()
axis(2, at=seq(-0.5,0.5,by=0.1), labels=c(-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1))
axis(1, at=seq(0,10), las=2, labels=labels)
segments(x.axis,m_police_BAGO$lower/2,x.axis,m_police_BAGO$upper/2)
points(x.axis,m_police_BAGO$coefficient/2, pch=21, bg=m_police_BAGO$pch)
abline(h=0, lty=3)
dev.off()

